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INTRODUCTION 


The  major  objective  of  this  article  is  to  examine  eigenvalue  equations 
expressed  in  terras  of  diraensionless  variables.  Our  application  is  to 
underwater  acoustic  ducts  for  which  the  square  of  the  index  of 
refraction  is  piecewise  linear.  For  this  case  the  eigenvalue  equation 
involves  the  Airy  functions  Ai ,  Bi,  Ai,  and  Bi  for  various  arguraents. 

We  consider  two  related  approaches  to  the  eigenvalue  problem.  The 
first  approach  is  the  usual  one,  while  the  second  approach  is  that  of 
this  article.  The  first  approach  is  to  determine  the  eigenvalues, 
i.e.,  the  mode  phase  velocity,  as  a  function  of  frequency  and  the 
profile  parameters.  Given  a  sound  speed  profile, one  iterates  the 
eigenvalue  equation  to  determine  the  mode  phase  velocity  as  a  function 
of  frequency.  This  approach  is  useful  when  treating  profiles  with  a 
large  number  of  layers. 

The  second  approach  is  to  express  the  eigenvalue  equation  in  terms  of 
dimensionless  variables.  These  variables  are  the  various  arguments  of 
the  Airy  functions  and  various  ratios  of  sound  speed  gradients.  The 
eigenvalue  solution  consists  of  a  set  of  these  dimensionless  variables, 
called  an  eigenvalue  set,  that  satisfy  the  eigenvalue  equation.  The 
frequency  and  phase  velocity  are  then  obtained  from  expressions  which 
involve  the  eigenvalue  set  of  diraensionless  variables  and  the  profile 
parameters.  The  solutions  are  valid  for  any  sound  speed  profile  that 
falls  within  the  layer  configurations  for  which  the  eigenvalue  equation 
is  formulated. 

We  note  that  the  two  approaches  use  exactly  the  same  eigenvalue 
equation.  In  the  first  approach  we  select  the  problem  variables 
(oroflle  parameters,  frequency,  phase  velocity)  and  convert  them  to 
dimensionless  mathematical  variables  which  are  then  used  in  the 
iterative  solution  of  the  eigenvalue  equation.  The  second  approach 


I 


reverses  these  steps.  We  first  determine  the  iterative  solution  of  the 
eigenvalue  equation  for  the  mathematical  variables.  These  mathematical 
variables  are  then  associated  with  the  physical  variables  for  the 
particular  profile  of  interest. 

As  we  shall  demonstrate,  the  second  approach  is  advantageous  in  the 
treatment  of  relatively  simple  duct  configurations  that  involve  only  a 
few  mathematical  variables.  Examples  are  given  that  illustrate  cases 
of  from  one  to  five  variables. 

The  second  approach  is  not  new.  For  example  it  is  used  in  Ref.  1,  which 
is  based  on  much  earlier  work  by  investigators  in  electromagnetic  propa¬ 
gation.  Reference  1  treats  a  bi-linear  surface  duct  for  which  there 

are  three  dimensionless  variables,  designated  as  Mx  ,  M  and  p. 

n 

Our  renewed  interest  in  the  second  approach  stems  from  the  studies  of 
Refs.  2  and  3.  In  Ref.  2  a  double  duct  with  nearly  coincident 
eigenvalues  was  constructed  and  investigated.  The  construction  was 
based  on  the  normal  mode  solutions  for  a  positive-gradient  surface  duct 
and  for  a  symmetric  refractive  duct.  For  these  simple  ducts  there  is 
only  one  dimensionless  variable  in  the  appropriate  eigenvalue  equation 
and  it  was  demonstrated  that  here  the  phase  velocity  could  be 
explicitly  represented  as  a  function  of  the  frequency,  profile 
parameters,  and  the  roots  of  Ai  and  Ai .  Furthermore  it  was  shown  that 
for  these  simple  ducts  the  phase  integral  results  of  ray  theory  could 
be  brought  into  congruence  with  mode  theory  by  the  use  of  non-integral 
values  of  mode  number. 

This  result  was  carried  forward  in  Ref.  3,  which  recommends  a  procedure 
for  testing  various  ray  theories  by  a  comparison  of  phase-integral 
results  with  the  exact  solutions  of  normal-mode  theory.  This  procedure 
pr'''iosed  the  use  of  non-integral  values  of  mode  number  such  that  the 
the  ray-  and  mode-theory  results  agree  in  the  high  frequency  limit. 


Section  I  presents  general  expressions  for  the  eigenvalue  formulation. 
Section  II  deals  with  the  unbounded  asymmetric  refractive  duct  for 
which  there  are  two  dimensionless  variables.  An  expression  for 
non-integral  mode  number,  which  brings  ray  and  mode  theory  into 
congruence,  is  developed.  Section  III  treats  an  asymmetric  refractive 
duct  with  a  surface  boundary.  For  this  case  there  are  three 
dimensionless  variables.  Section  IV  deals  with  a  surface  duct 
overlaying  a  refractive  duct.  For  this  general  case  there  are  five 
dimensionless  variables.  However  for  the  profile  configuration  of  Ref. 
2,  this  reduces  to  three  dimensionless  variables.  Section  V  outlines 
areas  for  further  investigation.  A  summary  is  provided  in  Section  VI. 

I.  GENERAL  EXPRESSIONS 

This  section  introduces  the  general  expressions  which  are  necessary  or 
useful  for  the  analysis  of  eigenvalues. 

The  sound  speed  in  each  layer  of  the  profile  is  expressed  as 

[c.j/c(z)]^  -  1  -  2y^(z  -  z^)/c^,  (1) 

where  c^,  z^,  and  are  the  sound  speed,  depth,  and  sound- 

speed  gradient,  respectively,  at  the  top  of  layer  i.  This  article  will 

only  treat  the  case  of  continuous  sound  speeds  at  layer  interfaces, 

i.e.,  the  sound  speed  at  the  top  of  layer  1  and  at  the  bottom  of  layer 

1-1  Is  the  same  with  a  value  of  c..  The  case  of  discontinuous  sound 

1 

speeds  is  tractable  but  leads  to  additional  complications. 

For  this  profile  the  unnormalized  depth  function  for  mode  n  and  layer  i 
may  be  written  as 

F  (Z)  =  A  Ai  (-f  )  +  8  Bi  (-C  .).  (2) 

m  m  m  m  m 
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Here  A  .  and  B  .  are  coefficients  which  are  independent  of  Z,  Ai 
m  m 

and  Bi  are  the  Airy  functions,  and  C.  is  given  by 


fni ( z) 


3  2  2  2  2 

a-j  (z  -  )  +  w  /c^  -  Xn  /^i  • 


(3) 


In  Eq.  (3) 


3  2  3 

a^  «  -2Y-iw  /Ci .  (4) 


The  quantity  x  is  known  by  several  names  i.e.,  the  mode  wave 
n 

number,  the  mode  eigenvalue,  and  the  separation  constant.  The  boundary 

conditions  and  the  interface  matching  conditions  form  a  system  of 

homogenous  linear  equations  in  the  coefficients  A^^  and  B^^.  The 

number  of  equations  is  equal  to  the  number  of  A^^  plus  the  number  of 

B  This  system  of  equations  has  a  non-trivial  solution  (non-zero) 
n  1 

if  and  only  if  the  determinant  of  the  coefficient  matrix  of  the  A  , 

ni 

and  B  ,  is  zero.  This  determinant  set  to  zero  is  the  eigenvalue 
ni 

equation.  The  x  are  the  values  of  mode  wave  number  for  which  the 
n 

determinant  is  zero. 


At  layer  interface  i  the  pressure  and  its  depth  derivative  must  be 

continuous.  These  matching  conditions  require  Airy  function  evaluation 

at  the  upper  interface  of  layer  i  and  the  lower  interface  of  layer 

i-1.  Thus  at  the  upper  interface  of  layer  i  the  Airy  functions  must  be 

evaluated  at  -x^  .  where 
i.i 


*1.1 


2  2 
Xn ) /a i . 


(5) 
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Equation  (5)  follows  immediately  from  Eq.  (3)  with  z  =  z^.  We  note 
that  Eq.  (3)  may  be  evaluated  in  terms  of  c  rather  than  z.  From  Eqs. 
(1)  and  (3)  it  follow,  that 


2  2 

Cni  ”  (w2/c2  -  Xf|)/a'f. 


(6) 


This  expression  yields  not  only  Eq.  (5)  but  also  the  lower  interface 
evaluation  of  layer  i-1  at  -x^  where 


-  2  2  2 

xi,i-l  =  («^/ci  -  ^n)/ai-l 


(7) 


At  layer  interface  i  the  continuity  of  pressure  leads  to 
An. i-1  Ai(-Xij.i)  +  Bpj.T  Bi(-Xi,i-1) 

-  Ani  Ai(-Xi  j)  -  Bni  Bi(-Xi^i)  -  0.  (8) 

The  continuity  of  depth  derivative  leads  to 

1  1 
An, i-1  ^i-1  Ai(-x.jj_-|)  +  Bnj_i  aj.i  Bi(-Xjj_i) 

We  now  examine  various  boundary  conditions.  Consider  first  cases  where 
the  layer  is  an  unbounded  half  space.  When  layer  i  is  an  unbounded 
half  space  with  positive  sound-speed  gradient,  the  solution  is 

i.e.,  the  coefficient  B  .is  zero.  Similarly  when  layer  i  is  an 
unbounded  half  space  with  negative  sound-speed  gradient,  the  solution  is 


F  Jz)  =  A  .  Ai 
nl  m 


s 


(11) 


In  this  latter  case  Eq.  (4)  cannot  be  used  directly  as  a^  is  evaluated 
in  terms  of  the  slope  and  sound  speed  at  the  upper  interface.  However, 
we  can  also  evaluate  a^  in  terms  of  the  slope  and  sound  speed  at  the 
lower  Interface.  From  Eq.  (1)  we  determine  that 

3  3 

dc/dz  =  c  .  (12) 

Thus  if  y.  is  the  slope  at  the  bottom  of  layer  i, 

10 


T 


io 


T/C, 


(13) 


Hence 


Ulo  “  ^'ui 


(14) 


evaluates  a^  in  terms  of  the  slope  and  sound  speed  at  the  lower 
Interface. 


Consider  next  the  case  where  interface  1  is  the  ocean  surface.  The 
condition  for  this  Interface  is 

%l'^>  -‘m  *'<-'l.l>  ^“al  61(-x,,,)  -0.  (15) 

Equations  (8)  to  (11)  and  (15)  will  allow  us  to  express  the  eigenvalue 
equation  for  the  duct  configurations  of  interest. 

Me  have  found  that  the  mode  phase  velocity  is  somewhat  easier  to 
Interpret  and  to  analyze  than  the  mode  wave  number.  The  mode  phase 
velocity  is  given  by 


c  =  (j/Rl  X. 
pn  n 


(16) 


The  modes  are  always  ordered  by  increasing  phase  velocity  i.e.,  mode  1 
has  the  lowest  phase  velocity,  mode  2  the  next  lowest  etc.  Equation 
(5)  may  be  expressed  in  terms  of  mode  phase  velocity  as 


M.i 


,  ,  ,  ,2  l/2,2/3  ,  ,-2/3 


Equation  (17)  may  be  solved  for  c^^  to  obtain 

,,  ,-2/3^, -1/2 

'pn  ■  'l<'  -f*>  ■ 


where 


2/3  2/3 

X  *  x^(-Y^  )  ”  "j  ^i 


(17) 


(18) 


(19) 


and 


-2/3  /->nv 

“j  *  •  "i.r 

The  simplest  configuration  is  the  single  positive-gradient  surface 
duct.  Here  Eqs.  (10)  and  (15)  both  apply  and  the  eigenvalue  equation 
reduces  to 

Ai(-x  )  =  0.  (21) 

'  f  ' 

Here  the  generic  x  represents  the  single  dimensionless  variable  of 

'  f  ^ 

this  profile  configuration.  The  eigenvalues  are  the  roots  of  the  Airy 
function  with  sign  changed.  Equations  (18)  to  (20)  with  i=l  then  gives 
an  explicit  expression  for  phase  velocity  in  terms  of  the  frequency, 
the  profile  parameters  c.^  and  ,  and  the  eigenvalues  (roots  of 
the  A1 ry  function) . 


II.  UNBOUNDED  REFRACTIVE  DUCT 


This  section  treats  the  solutions  for  a  two-layer  unbounded  refractive 
duct.  The  general  case  of  an  asymmetric  duct  is  treated  with  the 
symmetric  duct  as  a  special  sub  case.^ 

Figure  1  is  a  schematic  of  the  duct.  The  arrows  indicate  that  the  two 
layers  are  unbounded.  From  the  standpoint  of  the  normal  mode  solution, 
this  profile  can  be  characterized  by  three  parameters.  These  are  the 
axial  sound  speed  (c^),  the  gradient  at  the  axis  for  the  lower  layer 
(^.1).  and  for  the  upper  layer  (t^q)- 

For  the  upper  and  lower  layers  Eqs.  (10)  and  (11)  respectively  apply. 

At  the  layer  interface  (axis)  Eqs.  (8)  and  (9)  apply  to  yield 


*n0  ■  *r,l  ■  “• 


and 


AnO  ao  AK-xi^o)  '  ^nl  ^1  Ai(-xij)  =  0. 


We  now  let 


(22) 


(23) 


p  =  a^Z-a^,  (24) 

where  is  evaluated  by  Eq.  (14)  and  a.|  by  Eq.  (4).  These  expres¬ 
sions  lead  to 


p  = 


1/3 


(25) 
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We  now  let  x  =  The  eigenvalue  equation  may  then  be  expressed 

from  Eqs.  (22)  and  (23)  as 

I Ai(-p2x)  -  Ai(-x)  I 


1  1 
Ai(-p2x)  p  A1(-x) 


Expansion  of  Eq.  (26)  leads  to 


0. 


1 


(26) 


Gi(x,p)  =  p  Ai(-p2x)  Al(-x)  +  Ai(-x)  Ai(-p2x)  =  0.  (27) 

Thus  Eqs.  (18)  to  (20)  with  1»1  and  x^  ^  equal  to  the  root  of  Eq.  (27) 
yield  the  phase  velocity. 


Here  the  dimensionless  variables  are  p  and  x.  The  phase  velocity  Is 
obtained  from  Eq.  (18)  as  an  explicit  function  of  the  frequency,  the 
profile  parameters  c.^  and  y.|  .  and  the  dimensionless  eigenvalue  x. 

The  third  profile  parameter,  does  not  appear  explicitly  In 
Eqs.  (18)  to  (20).  However  It  appears  in  the  eigenvalue  equation 
through  the  variable  p  and  thus  Influences  the  value  of  x  which 
satisfies  Eq.  (27). 

We  note  that  the  solution  x  of  Eq.  (27)  applies  to  all  profiles  with 
the  given  value  of  p.  Moreover  one  eigenvalue  set  suffices  for  all 
frequencies.  The  advantages  of  this  characteristic  will  be  pointed  out 
later  in  this  section.  As  we  shall  see  later  for  more  complicated 
profile  configurations  each  eigenvalue  set  corresponds  to  a  single 
frequency. 

Consider  now  the  case  of  a  symmetric  duct  i.e.,  p=l .  Here  Eq.  (27) 
reduces  to 


1 

2  Ai(-x)  Ai(-x)  =  0. 


(28) 


Thus  the  ^  represent  the  roots  of  the  Airy  function  with  the  sign 
changed  and’the  roots  of  the  Airy  function  derivative  with  the  sign 
changed.  The  roots  of  the  derivative  correspond  to  the  odd  number  modes 
while  the  roots  of  the  function  correspond  to  the  even  number  modes. 

A  computer  routine  was  developed  to  solve  Eq.  (27)  by  Newton^s  method. 
The  procedure  starts  at  p=l  with  a  known  root  of  Ai(-x)  or  Ai(-x). 

The  value  of  p  is  decreased  by  successive  steps  of  ap  and  the 
solution  of  Eq.  (27)  obtained  for  each  step  by  the  iteration 


Xi+1  *  +  G-|  (x-j  ,p)/. 


aG-j 


ax 


(29) 


Xi 


where 
aG] 
ax 


Xi 


1  1 

(p3+l)  Ai(-Xi)  Ai(-p2xi)  -  px^  Ai(-x^)  Ai(-p2xi) 


(30) 


The  initial  estimate  of  x^  is  taken  to  be  the  solution  of  Eq.  (27) 
for  the  previous  value  of  p.  Once  the  iteration  of  Eq.  (29)  reaches 
the  desired  accuracy,  the  process  is  stopped,  p  is  decreased  by 
Ap  and  the  iteration  process  repeated. 


We  note  that  the  solutions  of  Eq.  (27)  for  0  <  p  <  1  suffice  for  all  p. 
If  the  result  of  Eq.  (25)  is  greater  than  1,  we  consider  the  reflection 
of  the  profile  about  the  axis.  For  this  configuration  p  <  1. 


Figure  2  presents  the  solution  of  Eq.  (27)  as  a  function  of  p  for  the 
four  smallest  roots  i.e.,  the  first  four  modes.  At  p=l  the  odd  roots 
correspond  to  the  negative  of  the  roots  of  the  Airy  derivative  and  even 
roots  to  the  negative  of  the  roots  of  the  Airy  function.  At  p=0,  the 


10 


roots  are  given  by  the  negative  of  the  roots  of  the  Airy  function. 

Thus  for  mode  1  the  values  of  x  increase  monotonical ly  from  for 

p«l  to  x-|  for  p*0.  For  modes  2  to  4  they  increase  respectively  from 
1 

x.|  to  x^,  from  x^  to  and  from  x^  to  x^.  When  p=0,  Eq.  (27) 
reduces  to 

1 

Ai(-x)  Ai(0)  =  0.  (31) 

Thus  the  results  of  Fig.  2  for  p*0  are  predictable. 

The  solutions  of  Eq.  (31)  are  the  same  as  Eq.  (21)  for  the  single 
positive-gradient  surface  duct.  There  are  two  distinct  configurations 
corresponding  to  p=0.  The  first  is  t.!*©  which  arises  from  an  isospeed 
half  space  below  the  axis.  The  second  is  y.jq=-*  which  arises  from 
the  limit  of  a  steep  negative  gradient  above  the  axis.  Although  the 
mathematical  eigenvalues  are  the  same  for  both  configurations,  the  phy¬ 
sical  results  are  quite  different.  For  y.|=0,  the  phase  velocity  of 
Eq.  (18)  reduces  to  the  axial  sound  speed.  For  Y.|Q=-*t  the  phase 
velocity  of  Eq.  (18)  depends  in  the  usual  manner  on  ,  f,  and  the 
roots  of  Ai(-x)»0. 

We  now  examine  approximate  solutions  to  Eq.  (27)  as  based  on  Taylor 
series  expansions  about  a  given  solution.  Let  x^  be  the  known  solu¬ 
tion  for  Pq.  Let 

x  =  Xq  +  ax  (32) 

be  the  solution  for  p  =  Pq  Ap.  Let 


2 

P  X 


2 

Pq  Xq  +  ay. 


(33) 


where 


2  2^ 

ay  =  ax  p  +  (p  -  pq)Xq.  (34) 

We  expand  Ai(-x)  and  a1(-x)  as  Taylor  series  in  ax  about  Xq,  and 

2^2  ^ 

Ai(-p  X)  and  Ai(-p  x)  as  Taylor  series  in  ay  about  p^  Xq. 

Consider  just  the  case  of  general  p^.  When  the  first  3  terms  of  the 
Taylor  expansions  in  ax  and  ay  are  substituted  in  Eq.  (27),  a  sur¬ 
feit  of  14  terms  result.  There  are  two  constant  terms  which  represent 
Eq.  (27),  evaluated  at  x^  and  p^,  and  which  sum  to  zero.  There 

are  two  terms  each  in  ax,  ay,  and  ax  ay.  There  are  three  terms 
2  2 

each  in  (ax)  and  (ay)  .  The  result  for  second-order  terms  is 
too  complicated  for  our  purpose  here.  If  we  retain  only  first  order 
terms  in  ax  and  ay,  we  obtain  as  a  correction 

2  2  2 

4x  =  (pq-p^)Xq  (p+Kpq)/1+p^+pI^(1+ppq) t  (35) 

where 

K  »  -  Xq  Ai(-pQXQ)  Ai(-XQ)/A](-pgXQ)  Ai(-Xq).  (36) 

Equations  (35)  and  (36)  then  give  the  value  of  x  for  p  in  a  neighbor¬ 
hood  about  Pq. 

However  If  the  expansion  point  is  taken  at  Pq=1  or  Pq=0,  the  expansion 

greatly  simplifies.  Consider  first  the  expansion  for  odd  modes  about 

111  1 
Pq=1 .  Here  =  x^ ,  where  x^  is  a  root  of  Ai(-x)=0.  Thus  both  Ai(-XQ) 

1  2 

and  ^'’("Pq^q)  3re  zero  and  the  14  terms  for  general  p^  collapse  to 

2 

only  4  terms  which  have  a  common  factor  of  .  Here  the 

ax  is  a  solution  to  the  quadratic 

2 

A(ax)  -t-  Bax  -*-0  =  0, 


(37) 


where 


A  ■  p(p2-p+l )/2, 

1  - 

B  *  X^p(p2-p+1 ) , 

and 

l2 

C  -  Xi  (p-1)  (Up2)/2. 

One  of  the  two  roots  of  Eq.  (37)  is 
is  obtained  by  choosing  the  sign  of 
as  positive.  This  choice  yields  Ax 
when  p*l . 


(38) 

(39) 

(40) 

spurious.  The  desired  root  of  Eq.  (37) 
2  1/2 

(B  -4AC)  in  the  quadratic  formula 
=0  as  the  correct  root  when  C=0,  i.e., 


In  the  case  of  even  modes  about  p«*1 ,  x-=x.,  where  x.  is  a  root  of 

0  0  1  1 
2 

Ai(-x)«0.  Here  both  Ai(-XQ)  and  Ai(-pQXQ)  are  zero  and  the  14  terms 
for  general  p^  collapse  to  only  2  terms  which  have  a  common  factor  of 
1  1  2 


A1  (-X 


0 


AX 


A1(-poXo).  Here  Ax  satisfies  a  linear  equation  which  leads  to 
2 

*  -  px^(p-l)/p  -  p  +  1.  (41) 


In  the  case  of  all  modes  about  Pq=0,  where  x^  is  a  root  of 

2 

A1(-x)*0.  Here  Ai(-XQ)»0.  However  Ai(-pQXQ)'=Ai(0) .  The  14  terms  for 

1 

general  pq  collapse  to  5  terms  which  have  a  common  factor  of  Ai(-XQ). 
Here  ax  is  a  solution  to  the  quadratic 

2 

A(ax)  +  Bax  +  C  «  0,  (42) 

where 

A  =  -Xi/2.  (43) 

B  =  X ( 1 +p3 )  ,  ( 44 ) 

C  =  p( 1 -p^Kx 1  )  , 


(45) 


and 


1 

K  -  A1(0)/Ai(0)  -  -0.7290112.  (46) 

The  desired  root  of  Eq.  (42)  is  again  obtained  by  choosing  the  sign  of 
2  1/2 

(B  -4AC)  in  the  quadratic  formula  as  positive.  This  choice  yields 
Ax*0  as  the  correct  root  when  C=0,  i.e.,  when  p=0. 

The  results  of  the  various  approximations  are  compared  with  the  exact 
solution  in  Tables  1  and  2  for  modes  1  and  2  respectively.  Column  2 
gives  the  exact  solution  of  Eq.  (27)  as  obtained  by  iteration.  Column  3 
presents  the  application  of  Eq.  (35)  for  the  expansion  point  Pq«0.5. 

The  value  of  Eq.  (37)  is  3.4268069  and  57.9192390  for  mode  1  and  2 
respectively.  Column  4  of  Tables  1  and  2  present  the  applications  of 
Eq.  (37)  and  (41)  respectively.  Column  5  presents  the  application  of 
Eq.  (42). 

Equation  (37)  remains  accurate  to  three  significant  digits  at  p*.65 
and  1s  the  most  accurate  approximation.  Equation  (41)  remains  accurate 
to  three  significant  digits  at  p«.85  and  is  the  second  most  accurate 
approximation.  In  contrast  the  accuracies  of  Eq.  (35)  and  Eq.  (42)  are 
somewhat  disappointing.  The  iterative  solution  of  Eq.  (27)  is  the 
method  of  choice.  If  accurate  algebraic  approximates  are  desired,  we 
recommend  that  the  quadratic  counterpart  of  Eq.  (35)  be  developed  and 
its  accuracy  about  various  expansion  points  assessed.  We  believe  that 
this  quadratic  counterpart  at  various  expansion  points  together  with 
Eqs.  (37),  (41)  and  (42)  can  provide  accurate  algebraic  approximations. 
This  approach  is  laborious  to  develop  but  is  straightforward. 

We  now  compare  the  solution  of  Eq.  (37)  with  that  of  Eq.  (41)  for 
values  of  p  near  1 .  Let 


p 


1  -Ap  . 


(47) 


We  expand  the  solutions  as  power  series  in  Lp.  The  general  form  is 

2 

ax  =  x(l)  Ap  [1  +  OAp  +  E(Ap)  ],  (48) 

where  x(l)  represents  the  solution  for  p»l .  In  the  case  of  Eq.  37, 

0-1/2  and  E— 1/2.  In  the  case  of  Eq.  (41),  0-0  and  E— 1 . 

We  now  assess  the  modified  phase  integral  result  of  Ref.  2.  If  we  use 
the  modified  n  of  Eqs.  (38)  and  (39)  of  Ref.  2  and  expand  Eq.  (6)  of 
Ref.  2  as  a  power  series  in  Ap  we  find  Eq.  (48)  holds  with  0-1/4 
and  E— 13/6.  Thus  we  see  that  the  modified  phase  integral  result 

agrees  with  the  result  of  Eqs.  (37)  to  (41)  to  first  order  in  Ap. 

Moreover  the  modified  phase  integral  result  lies  exactly  midway  between 
that  for  Eqs.  (37)  and  (41).  This  result  is  gratifying,  because  the 
modified  phase  integral  result  must  be  applied  to  both  odd  and  even 
modes.  The  result  properly  makes  a  compromise  by  "splitting  the 
difference"  between  the  approximation  for  odd  and  even  modes.  Our 
conclusion  is  that  for  value  of  p  near  1  the  modified  phase  integral 
result  of  Ref.  2  represents  a  fair  approximation  which  is  not  quite  as 
good  as  that  of  Eqs.  (37)  and  (41). 

We  now  present  a  modified  phase  integral  approach  which  results  in 
exact  values  of  phase  velocity  for  an  asymmetric  refractive  duct.  We 
set  Eq.  (19)  equal  to  Eq.  (6)  of  Ref.  2  and  solve  for  n.  The  result  is 

n  -  2x^^^  (l+p^)/3ir  +  1/2,  (49) 

where  x  is  the  solution  of  Eq.  (27)  for  the  given  p.  Note  that  when 
p-1 ,  Eq.  (49)  is  Identical  to  Eq.  (38)  or  (39)  of  Ref.  2.  When  p=0, 

Eq.  (49)  is  Identical  to  Eq.  (38)  of  Ref.  2.  Thus  Eq.  (49)  takes  on  the 
D'-oper  values  for  p=1  and  p=0. 


The  validity  of  £q.  (49)  was  further  checked  by  applying  the  method  to 
the  two  single  unbounded  ducts  treated  in  Ref.  4.  Figure  3  is  a  copy 
of  Fig.  31  of  Ref.  4.  The  circles  are  the  normal  mode  phase  velocities 
as  determined  for  each  duct  by  the  first  approach  described  here  in  the 
introduction.  The  dashed  and  solid  curves  represent  the  phase  integral 
result  for  the  upper  and  lower  ducts  respectively.  As  can  be  seen  there 
are  systematic  differences  between  the  ray  and  mode  results.  These 
differences  were  attributed  to  duct  asymmetry. 

To  apply  Eq.  (49)  we  first  evaluated  the  value  of  p  for  the  ducts. 

These  values  of  p  were  0.593636  for  the  upper  duct  and  0.806781  for 
the  lower  duct.  Table  3  presents  a  summary  of  results.  Column  1  is 
the  mode  number.  Column  2  is  the  modified  n,  which  applies  for  sym¬ 
metric  ducts,  as  determined  by  Eqs.  (38)  and  (39)  of  Ref.  2.  The 
entries  of  columns  3  and  5  are  the  roots  of  Eq.  (27)  for  the  upper  and 

lower  ducts  respectively.  Columns  4  and  6  are  the  corresponding  value 
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of  Eq.  (49)  as  determined  by  x  and  p  . 

The  curves  of  Fig.  4  are  the  counterparts  of  those  of  Fig.  3,  where  the 
non-integral  mode  numbers  of  Table  3  are  used  rather  than  the  integral 
values  of  the  phase  Integral  curves  of  Fig.  3.  We  have  overlaid  these 
curves  on  the  mode  data  of  Fig.  3  and  found  that  they  go  through  the 
centers  of  the  circles  within  the  plotting  accuracy.  Thus,  this  check 
demonstrates  that  for  the  asymmetric  refractive  duct  the  use  of  Eq. 

(49)  brings  the  results  of  the  phase-integral  method  of  ray  theory  into 
congruence  with  the  exact  normal  mode  solution. 

Figure  4  also  illustrates  an  advantage  of  the  dimensionless  variable  ap¬ 
proach.  Each  curve  of  Fig.  4  requires  only  one  set  of  iterations  of  the 
eigenvalue  equation,  i.e.,  Eq.  (27),  to  determine  the  x  for  the  given  p 
and  desired  mode  number.  In  contrast  the  circles  of  Fig.  3  were  ob¬ 
tained  by  the  usual  approach  in  which  the  phase  velocity  is  determined 
by  Iterating  the  eigenvalue  equation  for  each  desi'^ed  f'"eauehcv.  Thus 
the  generation  of  circles  ot  Fig.  3  requires  about  >00  t'mes  the  number 
of  sets  of  iterat'ons  as  did  the  generation  ot  the  rjrves 

I  f 


Observe  1n  F1g.  3,  the  larger  displacement  between  ray  and  mode  theory 

for  the  upper  duct  as  compared  to  the  lower  duct.  In  Ref,  4  we 

attributed  this  to  the  fact  that  the  p  for  the  upper  duct  was  smaller 

than  that  of  the  lower  duct.  We  now  know  that  this  conclusion  was  in 

error.  Consider  for  example  entries  4  and  6  of  Table  3  for  mode  1.  We 

see  that  n=l  (the  value  for  the  phase-integral  curves  lU  and  IL  of  Fig. 

3)  lies  closer  to  the  upper  duct  entry  than  to  the  lower  duct  entry. 

Why  then  does  the  phase-integral  curve  for  the  upper  duct  lie  further 

away  from  the  mode  theory  solution?  The  answer  lies  in  the  fact  that 
2/3 

(Y.j)  appears  as  a  factor  in  Eq.  (19).  This  factor  is  about 
4.3  larger  for  the  upper  duct  than  for  the  lower  duct.  Thus  if  the 
error  in  x  (n)  were  the  same  for  both  ducts  the  error  in  the  x  of  Eq. 
(19)  would  be  4.3  times  larger  for  the  upper  duct.  We  see  then  that 
the  larger  discrepancy  for  the  upper  duct  is  associated  with  a  larger 
gradient  rather  than  larger  p. 

We  close  this  section  with  the  presentation  of  Fig.  5  which  is  based  on 
the  double  duct  of  Ref.  4.  The  circles  represent  the  normal  mode 
results  for  modes  2  and  3  of  the  double  duct.  The  curves  are  the 
modified  phase  integral  results,  based  on  Eq.  (49),  for  mode  1  of  the 
upper  duct  and  mode  2  of  the  lower  duct. 

The  curves  are  good  approximations  to  the  double-duct  results  with  the 
exception  of  the  region  where  the  curves  cross.  Thus  we  see  that  the 
solutions  for  the  single  refractive  profile  of  Fig.  1  are  useful  in  the 
wider  context  of  double  ducts  with  a  surface  boundary.  The  same  conclu¬ 
sion  was  reached  in  Ref.  2  for  the  case  of  a  surface  duct  overlaying  a 
symmetric  refractive  duct.  However  this  conclusion  is  based  on 
numerical  examples.  One  of  the  remaining  tasks  of  this  article  is  to 
demonstrate  analytically  the  relationship  between  a  simple  unbounded 
duct  or  a  single-layer  surface  duct  and  more  complicated  profiles  such 
as  the  double  duct  of  Ref.  2. 


III.  REFRACTIVE  DUCT  WITH  SURFACE 


The  simplest  extension  of  an  unbounded  refractive  duct  is  to  bound  the  upper 
layer  with  a  free  surface.  Figure  6  presents  a  schematic  of  the  profile, 
which  is  characterized  by  four  parameters.  These  are  the  axial  sound  speed 
(c^),  the  gradients  at  the  axis  for  the  lower  layer  (y^).  *^*<1 
upper  layer  (y2o^’  sound  speed  at  the  surface  (c.|). 


We  let 


-  2  2  2 
X2,2  •  X  *  (‘^/C2~’^n)/®2- 

,22  2 
XI  1  «  y  -  (w2/ci-Xp)p2/a2, 


where 


p  =  ('^2'^“^20^ 


The  eigenvalue  matrix  may  be  expressed  as 


Ai(-y) 

Ai(-p2x) 


Bi(-y) 

Bi(-p2x) 


-Ai(-x)  -  0. 


Ai(-p2x) 


Bi(-p2x) 


pAi(-x) 


This  may  be  written  as 

G{x.y.p)  =  A1(-y)G2  B1(-y)G-|  =  0 

where  Gi  is  given  by  Eq .  (27),  and 


1  1 

G2  =  pBi(-D^x)  Ai(-x)  Ai(-x)  Bi(-p2x). 


(56) 


Equation  (55)  may  be  considered  as  expressing  the  eigenvalue  x  for  given 

values  of  the  parameters  p  and  y.  The  parameter  p  is  determined  by 

the  sound  speed  profile.  Although  we  are  free  to  choose  the  value  of  y, 

this  choice  specifies  the  frequency.  This  characteristic  may  be  demon- 

2 

strated  by  eliminating  from  Eqs.  (50)  and  (51)  and  solving  for 
frequency  to  obtain 

f  -  (p  x-y)  [l-(C2/c.|)  ]  ('Y2o)»  •  (57) 

The  steps  in  the  solution  are  as  follows: 

1.  We  choose  a  value  of  y. 

2.  The  eigenvalue  equation.  Eg.  (55),  is  solved  for  x,  for  the 
chosen  value  of  y  and  the  value  of  p  for  the  desired  profile. 

3.  The  frequency  is  determined  from  Eq.  (57),  with  the  use  of  y, 

X,  and  other  parameters  for  the  desired  profile. 

4.  The  phase  velocity  is  determined  from  Eq.  (18)  with  the  use  of 
X,  the  frequency  of  Eq.  (57),  and  the  parameters  c^  and 

i'or  the  desired  profile. 

If  we  have  available  the  solution  of  Eq.  (55)  for  all  values  of  p  and 
y  we  have  the  eigenvalues  for  all  frequencies  and  for  all  profiles  of 
the  generic  form  of  Fig.  6.  The  solution  of  Eq.  (55)  may  be  obtained 
by  the  iteration 

x^^l  =  Xt  +  [G/OG/ax)]  ,  (58) 


where 

aG/ax  =  Ai(-y)  (aG^/ax)  +■  Bi(-y)  (dG^/ax) 

Here  aG^/ax  is  given  by  Eq.  (30)  and 

1  1 

aG^/aG  ^  -  f  )  [Ai(-x)  Bi(-p2x)  -  px  Ai(-x)  6i(-p^x)j 


(59) 


(60) 


Figure  7  presents  the  first  five  roots  of  Eq.  (55)  as  a  function  of  y 

for  the  case  of  p=1 .  Consider  now  the  physical  interpretation  of  y. 

We  see  from  Eqs.  (51)  and  (16)  that  y=0  corresponds  to  c  =C-.  This 

P  1 

corresponds  to  the  case  of  a  ray  grazing  the  ocean  surface.  Positive 
values  of  y  represent  modes  with  phase  velocity  greater  than  the  surface 
i.e.,  rays  reflecting  from  the  ocean  surface.  Negative  values  of  y 
correspond  to  propagation  in  the  refractive  duct  with  no  reflection 
from  the  surface. 

Consider  the  case  of  large  negative  y.  From  Eq.  (57)  this  corresponds 
to  high  frequencies.  This  then  represents  the  case  of  strongly  trapped 
modes.  The  appropriate  asymptotic  expressions  for  Ai(-y)  and  Bi(-y) 
are  given  by  Eqs.  10.4.59  and  10.4.63  of  Ref.  5.  These  expressions 
lead  to 

Ai(-y)/Bi(-y)  ~  2  ^  exp  (-20,  (61) 

where 

3/2 

C  -  2(-y)  V3.  (62) 

Thus  Eq.  (61)  approaches  zero  for  large  negative  y  and  Eq.  (55)  reduces 
to  Eq.  (27).  Thus  we  have  reached  our  initial  goal.  For  high  frequen¬ 
cies  (strongly  trapped  modes)  the  eigenvalues  for  the  bounded  duct  of 
Fig.  6  approach  those  of  the  unoounded  duct  of  Fig.  1. 

Consider  now  the  solid  horizontal  lines  of  Fig.  7.  These  represent  the 
roots  of  Eq.  (28),  i.e.,  roots  of  the  Airy  function  and  its  derivative. 
Equation  (28)  is  the  solution  of  Eq.  (27)  for  the  special  case  of  p=l 
i.e.,  the  condition  of  Fig.  7.  We  see  that  the  roots  of  Eq.  (55) 
rather  rapidly  approach  those  for  the  unbounded  duct  as  y  is  decreased 
below  zero. 


2(1 


Consider  next  the  solid  vertical  lines  of  Fig.  7.  These  lines  are 
located  at  the  roots  of  A1(-y)=0.  We  see  that  Eq.  (55)  again  reduces 
to  Eq.  (27)  in  general  and  to  Eq.  (28)  for  p*l .  Thus  in  Fig.  7  each 
curve  crosses  the  vertical  solid  lines  at  the  horizontal  lines  which 
are  the  roots  of  Eq.  (28). 

Moreover  we  note  that  Eq.  (55)  reduces  to 

6^  =  0.  (63) 

when  B1(-y)=0.  For  the  special  case  of  p=l ,  Eq.  (63)  reduces  to 

1  1 

Bi(-x)  Ai(-x)  +  Ai(-x)  B1(-x)  =  0.  (64) 

With  the  use  of  the  Wronskian  for  the  Airy  functions  Eq.  (64)  may  be 
simplified  to 

Al(-x)  Bi(-x)  +  (2v)‘^  =  0.  (65) 

Equation  (65)  has  been  solved  by  iteration.  The  first  five  roots  are 
given  in  column  2  of  Table  4.  These  roots  are  plotted  as  dashed 
horizontal  lines  in  Fig.  7.  The  vertical  dashed  lines  correspond  to 
the  roots  of  Bi(-y)=0.  We  see  that  the  curves  of  Fig.  7  cross  the 
intersections  of  the  dashed  lines. 

We  now  see  the  advantage  of  the  dimensionless  variable  approach. 

Figure  7  provides  the  eigenvalue  for  the  first  five  modes  for  all 
profiles  of  the  generic  form  of  Fig.  6  with  the  restriction  that  the 
refractive  duct  is  symmetric. 

We  next  examine  Eq.  (57).  We  see  that  when 
2 

y  =  /=  X. 


(66) 


the  frequency  is  zero.  Figure  8  again  presents  the  eigenvalue  curves. 
The  slant  line  is  the  locus,  y=x,  i.e.,  Eq.  (66)  for  p=l .  Eigenvalues 
on  the  slant  line  represent  zero  frequency.  Eigenvalues  to  the  left  of 
the  slant  line  represent  real  frequencies  while  those  to  the  right 
represent  pure  imaginary  frequencies. 

The  points  of  intersection  of  Eq.  (66)  with  the  eigenvalue  curves  can 
be  expressed  analytically.  Substitution  of  Eq.  (66)  into  Eq.  (55) 
leads  to 

1  1 

2p  Ai(-p2x)  Bi(-p2x)  Ai(-x)  +  Ai(-p2x)  Ai(-x)  Bi(p-x) 

1 

+  Bi(-p2x)  Ai(-x)  Ai(-p2x)  =  0.  (67) 

For  the  special  case  of  p=l ,  Eq.  (67)  reduces  to 

1  1 

Ai(-x)  [3  Bi(-x)  Ai(-x)  +  Ai(-x)  Bi(-x)]  »  0.  (68) 

With  the  use  of  the  Airy  function  Wronskian  the  bracketed  expression  of 
Eq.  (68)  may  be  further  simplified  to 

1 

Ai(-x)  Bi(-x)  +  (Air)"^  =  0.  (69) 

This  equation  is  similar  to  Eq.  (65)  and  has  been  solved  by  iteration. 
The  solution  to  Eq.  (68)  then  consists  of  the  roots  of  Eq.  (69)  plus 
the  roots  of  Ai(-x)=0. 

The  first  five  roots  of  Eq.  (68)  are  given  in  column  3  of  Table  4. 

Roots  1,  2,  4,  and  5  are  the  first  four  roots  of  Eq.  (69),  whereas  root 
3  (and  6)  are  roots  of  Ai(-x)=0.  Thebe  roots  are  plotted  as  solid 
horizontal  lines  in  Fig.  8.  We  see  that  the  curves  of  Fig.  8  do  indeed 
cross  at  the  intersection  of  the  slant  line  with  the  horizontal  lines. 


Although  F1g.  7  only  presents  positive  values  of  x,  there  are  negative 
roots.  Figure  9  presents  the  continuation  to  negative  values  of  the 
curve  of  Fig.  7  for  the  lowest  order  root.  The  curve  crosses  zero  at 
the  first  root  of  Bi(-y)=0  as  previously  discussed  and  approaches  the 
first  root  of  Ai(-y)=0  as  a  vertical  asymptote.  This  latter 
characteristic  is  a  consequence  of  the  fact  that  both  terms  of  Eq.  (28) 
go  to  zero  as  x-»-«.  In  general  the  curve  of  order  n  will  cross 
zero  at  the  n'th  root  of  Bi(-y)=0  and  will  approach  the  n'th  root  of 
Ai(-y)=0  as  a  vertical  asymptote. 

We  will  not  present  the  counterpart  of  Fig.  7  for  other  values  of  p. 
However  each  value  of  p  has  a  horizontal  and  vertical  grid  of  lines 
which  bears  the  same  relationship  to  the  curves  as  shown  in  Fig.  7. 

The  vertical  grid  is  always  the  same  because  it  does  not  depend  on 
p.  The  solid  lines  of  the  horizontal  grid  are  the  solutions  Eq. 

(27).  For  example,  for  p=0.5  from  Table  1  and  2  we  see  that  the 
first  and  second  solid  horizontal  lines  would  lie  at  1.64  and  3.35 
respectively.  Thus  for  example  the  second  mode  curve  for  p*0.5  is 
asymptotic  to  3.35,  crosses  1.64  at  y=4.09,  and  approaches  a  vertical 
asymptote  at  y=5.52  for  x=-®.  This  last  property  holds  because  Eq. 

(27)  is  satisfied  for  all  p  with  the  possible  exception  of 

As  p  varies  from  1  to  0  the  solid  horizontal  line  for  root  n  of  Fig. 

7  moves  monotonical ly  up  from  its  present  position  to  the  n'th  root  of 
Ai(-x)=0.  This  is  readily  inferred  from  the  results  of  Fig.  2. 

In  the  case  of  Fig.  1  we  did  not  have  to  deal  with  p>l ;  whereas  we 
must  do  so  for  the  case  of  Fig.  6.  One  may  readily  verify  that  if  x  is 
a  root  of  Eq.  (27)  for  p,  then 


=  p 


2 


X 


X 
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is  a  root  for  p~^ .  For  example  the  root  of  Eq,  (27)  for  p=2  for  the 
first  order  root  is  x=0.25(l .64)=.41 .  Thus  as  p  varies  from  1  to  « 
the  solid  horizontal  lines  of  Fig,  7  move  monotonical ly  down  from  their 
present  positions  and  approach  x=0  as  a  limit. 

For  other  values  of  p  the  dashed  horizontal  lines  are  the  roots  of 
Eg.  (56).  We  note  that  the  solid  and  dashed  horizontal  lines  of  Fig.  7 
are  interleaved.  Thus  as  the  solid  lines  move  upward  with  decreasing 
p  so  will  the  dashed  lines.  Similarly  the  dashed  lines  will  move 
down  with  increasing  p.  There  may  be  some  question  regarding  the 
dashed  horizontal  line  at  zero.  The  elements  of  Eq.  (56)  were  expanded 
in  Taylor  series  about  the  point  x=0.  It  was  found  that  this  root 
moved  upward  for  p<l  and  downward  for  p>l .  Thus  for  p<l  or 

p>l  X  will  turn  negative  for  y  larger  or  smaller  than  the 

appropriate  root  of  Bi(-y)=0. 

For  other  values  of  p  the  horizontal  lines  of  Fig.  8  are  the  roots  of 

Eq.  (67),  The  x  curves  will  intersect  these  lines  along  the  line 
2 

y*p  X. 

In  summarizing  the  results  of  Fig.  7,  we  note  that  the  curves  are 
surprisingly  simple  and  well  behaved.  The  set  of  vertical  and 
horizontal  lines  provide  a  framework,  which  makes  their  characteristics 
quite  predictable.  The  dependence  of  frequency  in  Eq.  (57)  on  the 

eigenvalues  x  and  y  makes  the  problem  more  complicated  than  that  of  the 

unbounded  duct.  However  given  an  x  versus  y  eigenvalue  curve  we  can 
readily  generate  from  Eq.  (57)  the  frequency  associated  with  the 
desired  profile. 

Although  plots  giving  x  as  a  function  of  y  for  various  values  of  p 
will  provide  all  eigenvalues  for  the  profile  configuration  of  Fig.  6, 
we  decided  to  prepare  a  plot  giving  x  as  a  function  of  p  for  a  fixed 
value  of  y.  We  chose  as  our  example  y=0,  because  this  corresponds  to 
the  interesting  case  of  the  ray  which  grazes  the  ocean  surface  i.e., 
c  =c, .  Figure  10  presents  the  first  four  roots  of  Eq.  (55)  as  a 

p  1 


24 


function  of  p  for  the  case  of  y=0.  The  roots  were  again  obtained  by 
the  method  of  Eqs.  (50)  to  (60).  The  process  was  started  at  p*l 
using  the  values  of  x  from  Fig.  7  that  correspond  to  y=0.  The  curves 
of  Fig.  10  were  then  generated  by  moving  p  from  1  up  to  4  and  from  1 
down  to  zero. 


Our  first  reaction  to  F1g.  10  was  that  the  results  were  in  error.  We 

1 

note  that  these  results  go  to  the  roots  of  Ai(-x)=0  for  p=0,  whereas 
we  expected  them  to  go  to  the  roots  of  Ai(-x)=0  for  p*0.  If  we  set 
p=0  in  Eqs.  (27)  and  (50),  Eq.  (55)  can  be  simplified  to 


Ai(-x) 


Ai(-y)  Bi(0)  +  Bi(-y)  A1(0) 


=  0. 


(71) 


Thus  it  would  appear  that  Ai(-x)=0  provides  the  roots  for  Eq.  (71)  and 
indeed  it  does  except  for  our  unfortunate  choice  of  y=0. 


When  y=0  the  term  in  the  brackets  of  Eq.  (71)  is  zero.  This  opens  the 

possibility  of  roots  other  than  Ai(-x)=0  for  p=0.  It  became  evident 

that  a  more  sophisticated  analysis  was  required.  We  let  p=c, 

2 

expressed  the  various  Airy  functions  with  argument  -p  x  as  two-term 
Taylor  series  expanded  about  zero,  substituted  these  series  in  Eq. 

(55),  and  collected  terms  of  various  order  in  c.  The  zero-order  term 
was  the  left  side  of  Eq.  (71),  which  is  zero  for  y=0  as  previously 

discussed.  The  first-order  term  was  2cAi(0)  Bi(0)  Ai(-x).  There  was 

3  ^  ^ 

no  second  order  term.  The  third-order  term  was  -2c  xAi(O)  Bi(0)  A1(-x). 

4  2  "I 

The  fourth-order  term  was  2c  x  Ai(0)  Bi(0)  Ai(-x).  We  see  that  Ai(-x)=0 

is  Indeed  the  limiting  root  as  p=c->0,  because  this  results  in  all 

3 

items  of  order  through  c  to  be  zero. 


In  o'-de’'  to  verify  the  analysis  of  Eq.  (71).  we  examined  the  first  root 
of  Ed.  (55)  as  a  function  of  p  for  the  case  of  y=-l .  Here  the  root 
went  to  2.338  for  p=0  as  predicted  by  Eq.  (71).  Further  analysis  is 
beyond  the  scope  of  the  present  article.  It  would  be  of  interest  for 
examp’e  to  examine  the  behavior  for  fixed  values  of  y  near  zero  for 


values  of  p  near  zero.  Would  the  curves  for  y=+c  be  nested  about 
the  curves  of  Fig.  10  for  small  p?  If  so  how  do  they  approach  a 
different  limit  at  p=0? 

We  note  that  Fig.  10  exhibits  another  characteristic  that  we  discussed 
in  connection  with  Fig.  7,  i.e.,  x  increases  as  p  decreases  below  1 
and  decreases  as  p  increases  above  1. 

IV.  DOUBLE  DUCT  PROFILE 


This  section  addresses  one  general  form  of  double  duct  profile.  A 

brief  discussion  of  the  application  of  the  method  to  more  general 

profiles  is  presented.  Figure  11  presents  the  schematic  of  a  double 

duct  profile  which  consists  of  a  surface  duct  overlaying  a  refractive 

duct.  The  duct  can  be  characterized  by  six  parameters.  These  are  the 

surface  sound  speed  (c^),  the  barrier  sound  speed  (c^),  the  axial 

sound  speed  (c  )  and  the  gradients  at  the  top  of  the  three  layers 
0 

i.e.,  ,  y^,  and  Other  gradients  of  interest  such  as  y^q  and 

Y30  fwy  be  derived  from  the  given  parameters  with  the  use  of  Eq.  (13). 
The  surface  sound  speed  may  be  larger  or  smaller  than  the  axial  sound 
speed.  However  both  c  and  c  must  be  less  than  c  . 

I  0  C 


There  are  five  mathematical  variables  which  we  define  as  follows: 


,  2 

X3^3  =  X  =(«2/c3 

2  2 

-  w2/cp)/a3. 

(72) 

,  2 

^2,2  ~  ^  =(u)^/C2 

,  2  2 
-  w//Cp)/a2, 

(  73) 

2 

XI  j  =  y  =(^2/ci 

-  2  2 
-  u)2/Cp)/ai , 

(  74  ) 

1/3 

p  = 

"5  ; 

ana 


/■  _  -y  ' 
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The  eigenvalue  matrix  may  be  expressed  as 


A1(-y)  Bi(-y)  0 

0 

0 

2  2 

Al(-piw)  Bl(-piw)  -A1(-w) 

-B1(-w) 

0 

(77) 

1  2  1  2  T 

Al(-piw)  B1(-piw)  PlAI(-w) 

1 

PlBI(-w) 

0 

-  0. 

0  0  A1(-p2x) 

B1(-p2x) 

-A1(-x) 

0  0  a1(-p2x) 

b1(-p2x) 

paI(-x) 

This  may  be  written  as 

A1(-y)[G^W^-G2W2]  +  BK-ylCG^W^-G^W^] 

=  0. 

(78) 

where 

2  1  12 

Wt  =  ptBI(-piw)  B1(-w)  +  Bl(-piw)  B1(-w). 

(79) 

2  1  12 

Wj  =  p-|B1(-p-,w)  A1(-w)  +  B1(-piw)  A1(- 

-w) . 

(80) 

2  1  12 

W3  =  piAI(-piw)  B1(-w)  A1(-piw)  B1(-w), 

(81) 

2  1  12 

W4  -  ptAI(-piw)  A1(-w)  +  Al(-piw)  Ai(-w). 

(82) 

and  and  G^  are  given  by  Eqs.  (27)  and  (56)  respectively. 

In  carrying  forward  the  solution  for  the  profile  of  Fig.  ll,  1t  Is 
convenient  to  discuss  It,  as  well  as  the  simpler  profiles  already 
t'^eated,  1n  the  context  of  a  general  profile  consisting  of  n  interfaces 
and  m  ooundarles.  For  such  a  profile  there  are  ?n-HTH-i  profile 

parameters  and  the  two  physical  parameters  f  and  c^  for  a  total  of 
2n-HTH-3  i/ariaDles.  In  the  first  approach  2n-t-rrH-2  ya''iables  are 

considered  independent  while  the  dependent  variable,  c  .  is 

P 

constrained  by  the  eigenvalue  equation. 


In  the  second  approach  the  number  of  introduced  mathematical  variables 
is  2n+m.  This  is  the  number  of  interface  and  boundary  conditions  the 
eigenvalue  matrix  must  satisfy.  Moreover  there  are  introduced  2n+m 
constraints  which  define  the  mathematical  variables  in  terms  of  the 
profile  and  physical  parameters.  We  treat  2n-Hn-i  of  the  mathematical 
variables  as  independent  with  the  remaining  one  satisfying  the  eigen¬ 
value  equation.  Thus  of  the  2n+flH-3  profile  and  physical  variables, 

2n-Hii  are  dependent  while  the  remaining  three  are  independent. 

Consider  in  this  context  the  simpler  profiles,  already  discussed.  In 
the  case  of  the  single-layer  surface  duct  the  three  independent 
variables  are  c, ,  y, ,  and  f  while  the  dependent  variable  is  c  as  given 

1  1  p 

by  Eq.  (18).  In  the  case  of  the  profile  of  Fig.  1,  the  three  indepen¬ 
dent  variables  are  again  c^,  and  f.  The  two  dependent  variables  are 

y,^  as  constrained  by  Eq.  (24),  and  c  as  given  by  Eq.  (18).  In  the 

i  0  p 

case  of  the  profile  of  Fig.  6,  the  three  independent  variables  are  c^, 

^2’  ^20'  three  dependent  variables  are  y^  as  given 

by  Eq.  (52),  f  as  given  by  Eq.  (57),  and  c  as  given  by  Eq.  (18). 

P 

We  are  now  ready  to  proceed  with  the  case  of  Fig.  11  for  which  n=2  and 
m*l .  We  thus  must  select  three  Independent  and  five  dependent 
variables.  For  the  independent  variables  we  select  c  c  ,  and  y 
The  dependent  variable,  y  ,  is  constrained  by  Eq.  (75)  and  the  relation- 

3  ^ 

ship  "^30'  dependent  variable,  f,  is  constrained 

by 

f  =  (p  x-w;  [1-(C2/C2)  ] 

2 

Equation  (83)  was  obtained  by  eliminating  c  from  Eqs.  (72)  and  !'’3) 

p 

and  solving  for  f.  Equation  (83)  is  the  same  as  Eq.  \57j  with  the 
s  jt)St  ’  tut  i  on  0^  counterpart  parameters  bvven  ihis  ‘’'equency  ore  may 


determine  the  dependent  variable  from  Eq.  (18)  using  the  eigen¬ 
value  X  and  the  appropriate  parameters  associated  with  the  third  Inter¬ 
face  of  the  profile.  We  next  determine  the  dependent  variable,  c.| , 
from  the  constraint 

2/3  -2/3  -2/3  -2  -2 

Cl  *  (Y20  f  »  C2  y  +  Cp 

Equation  (84)  was  obtained  by  solving  Eq.  (74)  for  c.| .  Our  fifth  and 
final  dependent  variable  y.|  1s  constrained  from  Eq.  (76)  and  the 
relationship  -  (c^/c^I^Yjq- 

We  see  that  the  build  up  of  the  number  of  dependent  variables  with 
Increasingly  more  complicated  profiles  puts  severe  limitations  on  the 
second  approach.  In  fact  the  nature  of  Eq.  (84)  renders  the  second 
approach  essentially  useless.  The  basic  idea  behind  the  second 
approach  was  to  determine  eigenvalues  in  terms  of  the  mathematical 
variables  and  as  Independently  of  profile  parameters  as  possible.  The 
problem  with  the  profile  of  Fig.  11  Is  that  we  cannot  express  c.|  in 
terms  of  mathematical  variables  In  a  useful  manner.  For  example  a 
given  eigenvalue  set,  together  with  other  profile  parameters,  produces 
a  specific  value  of  c.^  which  In  general  will  not  be  useful.  On  the 
other  hand  we  cannot  readily  determine  other  variables,  which  will 
produce  a  desired  value  of  c^  as  these  variables  are  not  Independent 
of  each  other.  For  example  x  Is  a  function  of  y,  f  Is  a  function  of 

c.,  and  X.  and  c  Is  a  function  of  f. 

2  P 

At  this  stage  we  should  note  that  our  choice  of  the  three  independent 
variables  as  c^,  c^,  and  Yj  is  somewhat  arbitrary.  For  certain 
parametric  studies  of  the  theory  it  may  be  of  interest  to  consider 
other  sets  of  three  variables  as  independent.  However  the  use  of  other 
sets  will  not  solve  the  basic  problem  of  the  profile  of  Mg.  11.  The 
one  boundary  and  two  interfaces  lead  to  five  dependent  varidbies  which 
are  too  many  for  the  approach  to  cope  with 


The  problem  becomes  even  worse  for  example  if  the  profile  has  four 
distinct  sound  speeds  at  various  interfaces  or  boundaries.  Here  two  of 
these  sound  speeds  must  be  dependent  variables  and  will  be  subject  to 
various  constraints  such  as  Eq.  (84).  We  note  at  this  point  that  the 
problem  with  the  second  approach  stems  from  the  presence  of  more  than 
two  distinct  sound  speeds  at  interfaces  or  boundaries.  If  the  profile 
is  limited  to  two  distinct  sound  speeds  then  the  second  nnethod  is 
useful  and  is  similar  to  the  case  of  Sect.  Ill  where  the  x  and  y 
variables  are  associated  with  the  two  sound  speeds.  There  will  be  a 
suite  of  p  constraints  rather  than  the  single  p  of  Sect.  III. 

However  these  constraints  pose  no  fundamental  problem  to  the  second 
method,  other  than  to  complicate  the  eigenvalue  equation.  Examples  of 
this  will  be  given  later. 

Consider  now  the  case  of  strongly  trapped  modes  for  the  case  of  Fig. 

11.  From  Eq.  (83)  we  see  that  high  frequencies  correspond  to  large 
negative  w.  If  we  replace  each  Airy  function  in  Eqs.  (79)  to  (82)  by 
its  first  asymptotic  term  and  evaluate  for  large  negative  w  we  find  that 


^0 


For  large  negative  values  of  w,  Eq.  (90)  reduces  to 


Ai(-y)  =■  0.  (91) 

Now  is  Eq.  (27),  the  solution  for  the  unbounded  refractive  duct. 

A1(-y)»0  Is  the  solution  for  the  positive-gradient  surface  duct.  Thus 
at  high  frequencies  the  eigenvalues  for  the  profile  configuration  of 
Fig.  11  become  the  composite  of  those  for  the  surface  duct  and  the 
unbounded  refractive  duct. 

Figure  5  Illustrated  numerically  that  at  high  frequencies  the 
eigenvalues  for  a  double  duct  were  related  to  those  of  the  single 
ducts.  The  derivation  of  Eq.  (91)  has  demonstrated  this  result 
analytically.  The  result  is  not  new  but  the  manner  in  which  it  arises 
is  of ^interest.  The  result  arises  from  the  asymptotic  behavior  of  Ai , 
Bi,  Ai ,  and  Bi  and  the  particular  location  of  these  elements  in  the 
eigenvalue  matrix  of  Eq.  (77). 

We  have  already  demonstrated  that  the  second  approach  is  not  practical 
for  the  general  configuration  of  Fig,  11.  Consider  now  Table  5  which 
outlines  the  various  cases  that  arise  when  the  six  profile  parameters 
are  related  by  various  conditions.  Case  numbers  are  assigned  in  column 
1  for  ease  of  identification.  Column  2  gives  the  condition  and  column 
3  lists  the  mathematical  variables. 

Case  1  is  the  general  case  with  six  independent  profile  parameters  and 
five  mathematical  variables.  Case  2  to  4  represent  the  equality  of 
various  pairs  of  gradients.  Case  2  eliminates  p  and  case  3 
eliminates  p.^  .  The  condition  of  case  4  leads  to  the  constraint 

pPt  =  C  ^ /c^ .  (92) 

However  tn^s  does  not  eliminate  any  oi  the  mathemai  !ca  va'''at  es 


Case  5  1s  significant  because  It  leads  to  elimination  of  y  as  a 

variable  and  to  the  necessity  for  Eq.  (84)  as  a  constraint.  To 

demonstrate  this  we  note  that  p=a„/a„  and  p,=a^/a,.  Hence 

3  2  12  1 

Equation  (93)  is  true  for  the  general  profile.  However  if  c  =c  , 

I  J 

then  Eq.  (74)  may  be  written  as 

y  •  (pp^i^x  (94) 

and  y  is  eliminated  by  an  expression  involving  three  of  the  other 
mathematical  variables.  We  see  that  the  second  approach  now  becomes 
viable.  The  situation  is  akin  to  that  of  Sect.  Ill  except  we  must  deal 
with  two  gradient  ratios  rather  than  one. 

Case  6  is  the  combination  of  cases  2  and  5  and  eliminates  y  and  p. 

The  profile  investigated  in  Ref.  2  is  of  this  type.  Case  7  is  the 

combination  of  cases  3  and  5  and  eliminates  y  and  p^ .  Case  8  is 
the  combination  of  cases  4  and  5  and  eliminates  y  and  p^ .  Equation 
(93)  holds  and  reduces  to 

PP^  =  1  (95) 

for  c^=c^.  Thus  Eq.  (95)  may  be  used  to  eliminate  p^ . 

Case  9  is  the  combination  of  cases  2  anti  3  anti  eliminates  p  and 

p^  .  Cases  10  and  11  are  the  combination  of  case  4  with  cases  2  and 

3  respectively.  Neither  results  in  any  change  from  cases  2  and  3 
because  Eq.  (93)  does  not  reduce  the  variables  unless  c^=c^. 

Finally  case  12  is  the  combination  of  cases  2  to  5  and  eliminates  all 
Put  the  X  and  w  variables,  we  note  that  given  any  three  of  the  four 
conditions  case  12  implies  the  fourth  condition.  Thus  there  are  no 
d'Stinct  cases  involving  three  conditions. 


In  closing  this  section  we  note  that  the  vital  feature  of  profile 
simplification  in  Table  5  is  the  elimination  of  the  variable  y.  The 
apparent  elimination  of  gradient  ratios  by  making  them  equal  to  unity 
is  not  really  a  significant  simplification.  The  ratio  of  unity  still 
remains  as  a  constraint  on  the  profile  parameters  and  the  eigenvalue 
equation  is  not  really  less  difficult.  For  example  the  solution  of  Eq. 
(78)  is  not  significantly  more  difficult  for  p-0.5  than  it  is  for 
p»l .  The  elimination  of  gradient  ratios  may  actually  have  an  adverse 
effect.  For  example  we  suspect  that  case  12  of  Table  5  leads  to 
degenerate  eigenvalues,  because  of  the  symmetries  involved. 

V.  AREAS  FOR  FURTHER  INVESTIGATION 

Much  of  the  work  outlined  here  is  concerned  with  the  comparison  of 

various  ray  theory  approaches  with  normal  mode  theory  by  means  of  the 

phase  and  group  velocity.  This  is  the  major  thesis  of  Ref.  3,  which 

was  written  before  the  present  article  was  begun.  We  have  shown  that 

at  high  frequencies  the  mode  theory  phase  and  group  velocities  for 

various  bounded  ducts  goes  to  the  solution  for  simple  unbounded  ducts. 

Moreover  through  the  use  of  non-integral  n  we  can  make  ray  theory  phase 

and  group  velocities  agree  with  the  modal  result  for  the  simple 

unbounded  ducts.  The  idea  is  to  use  non-integral  n  in  the  phase 

Integral  method  which  would  utilize  “corrected"  expressions  for  R,  T, 

E.,  and  E  as  specified  by  the  ray  theory  approach  under  test.  The 
1  0 

use  of  non-integral  n  would  be  superior  to  the  use  of  integral  n 
because  the  ray  theory  result  would  now  go  to  the  exact  high  frequency 
limit. 

There  are  two  general  extensions  of  the  present  article  which  would  be 
useful  in  the  comparison  of  ray  theory  approaches.  We  have  presented 
the  solution  for  a  positive-gradient  duct  with  a  free-surface 
Doundary.  The  first  extension  is  to  determine  non-integral  values  of  n 
tnat  will  make  ray  theory  exact  for  a  negative-gradient  duct  with  a 


rigid-bottom  boundary.  In  this  case  the  non-integral  n's  are  a  function 
of  the  roots  of  Ai  rather  than  Ai .  A  second  extension  is  to  develop 
and  analyze  expressions  for  group  velocity  making  use  of  the 
dimensionless  mathematical  variables. 

Reference  3  recommended  three  ray  theories  for  comparison.  We  will 
only  outline  them  here.  There  are  two  important  updates  to  Ref.  3. 

The  first  is  that  we  now  have  available  Eq.  (49)  for  treating 
asymmetric  ducts.  The  second  is  that  we  recognize  that  the  second 
method  of  determining  eigenvalues  is  feasible  as  long  as  no  more  than 
two  distinct  sound  speeds  are  present  in  the  profile  model. 

One  of  the  ray  theories  cited  in  Ref.  3  is  associated  with  the 
introduction  of  medium  attenuation  through  the  use  of  complex 
coefficients  in  the  sound  speed  model.  We  first  note  that  the  solution 
of  Sect.  II  remains  valid  for  complex  parameters.  Thus  the  phase 
integral  method  with  non-integral  n  gives  the  identical  solution  for 
mode  attenuation  as  does  the  normal  mode  theory.  The  implementation  of 
this  approach  will  require  the  solution  of  Eq.  (27)  for  complex  p. 
Furthermore  the  n  of  Eq.  (49)  will  have  an  imaginary  component.  There 
are  two  other  ray  approaches  which  can  be  tested  against  this 
solution.  The  first  approach  integrates  the  local  attenuation  along 
the  ray  path  and  divides  by  R  to  obtain  the  attenuation  coefficient. 

The  second  approach  uses  rays  with  complex  phase  velocity.  Here  the 
attenuation  coefficient  is  given  by  wImT/R. 

A  second  ray  theory  for  testing  is  that  of  Ref.  5  which  treats  free  or 
rigid  boundaries.  The  profile  of  Fig.  6  can  be  used  to  test  a  free 
boundary.  The  left  profile  of  Fig.  12  can  be  used  to  test  both  free 
and  rigid  boundaries. 

A  third  -ay  theory  for  testing  is  that  of  Ref.  6  which  treats  rays 
turning  near  or  barely  penetrating  a  relative  maximum  in  sound  speed. 


Figure  11  with  c  =c  represents  a  semi -bounded  profile  suitable  for 

I  w 

such  testing.  The  right  profile  of  Fig.  12  presents  a  simpler  profile 
for  such  a  test. 

We  now  turn  to  further  work  not  concerned  with  testing  ray  theories. 

We  first  note  that  there  remain  some  interesting  questions  about  Sect. 
III.  What  is  the  significance  of  pure  imaginary  frequency?  What  is 
the  behavior  for  small  y  and  p  as  previously  discussed?  We  call 
attention  to  another  feature  of  Fig.  8.  Suppose  y=1.0.  Here  the 
smallest  eigenvalue  does  not  correspond  to  a  real  frequency.  Thus  the 
second  eigenvalue  must  correspond  to  the  first  mode.  However  as  y  is 
decreased  the  second  eigenvalue  must  correspond  to  the  second  mode  when 
the  frequency  of  the  first  eigenvalue  turns  real.  This  implies  that 
the  eigenfunction  for  the  second  mode  must  have  one  node  at  y=1.0  and 
this  must  change  to  two  nodes  as  y  is  decreased.  Consider  the  more 
extreme  case  of  the  fifth  eigenvalue  curve  in  Fig.  8.  At  y=-3,0  this 
is  the  fifth  mode  with  5  nodes  in  the  eigenfunction.  As  y  is  increased 
the  mode  number  (and  number  of  nodes)  will  decrease  by  one  each  time 
one  of  the  lower-order  eigenvalues  crosses  the  line  y=x. 

We  also  recommend  that  the  counterpart  of  Fig.  7  be  generated  for  the 
profile  of  Fig.  11  with  This  would  not  only  represent  a 

more  complicated  application  of  the  second  method,  but  would  shed  some 
analytic  Insight  into  some  of  the  unresolved  problems  of  Ref.  2. 

VI.  SUMMARY 

Two  general  approaches  to  eigenvalue  problems  have  been  considered.  In 
the  first  or  usual  approach  the  dimensionless  mathematical  variables 
are  evaluated  numerically  in  items  of  the  physical  and  profile 
parameters  The  eigenvalue  matrix  is  then  iterated  by  numerical 
methods  to  determine  the  mode  phase  velocity  in  terms  of  the  frequency 


and  profile  parameters.  In  the  second  approach,  which  is  the  main 
subject  of  this  article,  the  eigenvalue  equation  is  initially  solved  in 
terms  of  the  mathematical  variables  for  some  generic  profile 
configuration.  The  mode  phase  velocity  is  then  evaluated  in  terms  of 
the  frequency,  mathematical  eigenvalues,  and  the  parameters  for  any 
desired  profile  of  the  generic  configuration. 

The  simplest  generic  configurations  only  involve  profiles  with  two 
parameters.  There  are  two  configurations  of  this  type.  The  first  is  a 
surface  duct  consisting  of  a  positive-gradient  half  space  bounded  above 
by  a  free  surface.  The  second  is  a  negative-gradient  half  space 
bounded  below  by  a  rigid  bottom  surface.  Here  there  is  one  variable  in 
the  eigenvalue  equation  whose  solution  is  the  roots  of  Ai(-x)=0  for  the 
free  surface  and  Ai(-x)=0  for  the  rigid  surface.  The  phase  velocity 
may  then  be  solved  in  terms  of  these  Airy  function  roots  and  the 
independent  variables  of  frequency  and  two  profile  parameters. 

The  next  more  elaborate  configuration  involves  a  profile  with  three 
parameters.  This  is  an  unbounded  refractive  duct.  Here  there  are  two 
mathematical  variables  x  and  p  (a  ratio  of  gradients).  Here  the 
eigenvalue  x  is  solved  as  a  function  of  p.  These  eigenvalues  have 
been  computed  for  the  first  four  modes  for  0<p<l .  Values  for 
P>1  can  be  simply  expressed  in  terms  of  the  solution  of  p  \ 

Once  the  value  of  x  has  been  obtained  the  solution  proceeds  as  for  the 
simplest  configuration. 

The  next  more  elaborate  configuration  involves  two  different  sound 
speeds.  A  detailed  evaluation  has  been  carried  out  for  one  such 
configuration.  This  is  a  refractive  duct  bound  above  by  a  free  surface 
but  unbounded  below.  Here  the  profile  has  four  parameters  while  there 
are  three  mathematical  variables.  These  are  x  and  y,  associated  with 
tne  two  respective  sound  speeds,  and  the  ratio  of  axial  gradients  p. 
Here  the  eigenvalue  x  is  solved  as  a  function  of  y  for  a  fixed  value 

p . 


These  eigenvalues  have  been  computed  for  the  first  five  roots  for 
p=1 .  The  behavior  for  other  values  of  p  can  be  readily  inferred. 

For  any  configuration  with  two  distinct  sound  speeds  the  frequency  can 
no  longer  be  treated  as  an  independent  variable.  It  is  a  dependent 
variable  which  is  a  function  of  x  and  y  and  the  profile  parameters. 

Once  the  frequency  has  been  determined  as  a  dependent  variable  the 
solution  proceeds  as  previously  discussed. 

The  second  approach  cannot  be  carried  out  when  there  are  three  distinct 
sound  speeds  at  interfaces  or  boundaries  of  the  generic  profile 
configuration.  This  is  demonstrated  for  the  case  of  a  surface  duct 
overlaying  a  refractive  duct.  Here  the  three  sound  speeds  are  that  at 
the  surface,  at  the  duct  axis,  and  at  the  barrier  between  the  ducts. 

The  problem  here  is  that  all  three  sound  speeds  cannot  be  specified 
independently.  One  of  the  sound  speeds  must  be  expressed  as  a 
dependent  function  of  the  mathematical  eigenvalues  in  order  to  satisfy 
the  constraints  between  profile  parameters  and  mathematical  variables. 

However  the  second  approach  is  viable  for  degenerate  profile 
configurations  where  the  interface  and  boundary  sound  speeds  reduce  to 
only  two  distinct  values.  An  example  is  a  double-duct  configuration, 
just  referred  to,  for  which  the  surface  and  axial  sound  speeds  are  the 
same.  Other  examples  are  a  refractive  duct  bounded  above  and  below  by 
boundaries  at  the  same  sound  speed  or  a  positive-gradient  surface  duct 
overlaying  a  negative-gradient  bottom  reflected  duct  with  the  surface 
and  bottom  sound  speeds  the  same.  For  these  degenerate  configurations 
the  solution  proceeds  in  essentially  the  same  manner  as  described 
previously  for  the  refractive  duct  bound  above  by  a  free  surface  and 
unbounded  below.  The  chief  difference  is  that  the  eigenvalue  equation 
is  more  complicated  and  may  involve  several  ratios  of  gradients,  i.e., 
one  ratio  for  each  profile  interface 


The  second  approach  has  several  advantages  over  the  first  approach. 

One  advantage  is  that  the  second  approach  solves  the  eigenvalue 
equation  for  a  generic  profile  configuration.  These  eigenvalues  may 
then  be  used  to  determine  phase  velocity  for  any  desired  member  of  the 
profile  configuration.  In  contrast  the  first  approach  treats  each 
profile  and  frequency  as  a  distinct  problem  and  is  in  effect  solving 
the  same  problem  over  and  over  again.  Consider  for  example  the 
unbounded  refractive  duct.  The  second  approach  requires  one  set  of 
iterations  for  each  of  the  desired  values  of  p.  The  first  approach 
requires  a  set  of  iterations  for  each  frequency  desired  for  a  given 
profile  and  the  entire  process  must  be  repeated  for  each  profile 
desired. 

Another  advantage  of  the  second  approach  is  the  relative  ease  with 
which  the  eigenvalues  can  be  determined.  The  presented  numerical 
examples  exhibit  a  simple  behavior  with  a  wealth  of  mathematical 
properties  which  can  be  used  to  interpret  the  behavior  of  the 
solution.  An  excellent  example  is  the  refractive  duct  with  surface 
boundary.  Here  the  eigenvalue  curves  go  through  the  lattice  points  of 
straight  horizontal  and  vertical  lines  which  occur  at  various  roots  of 
the  Airy  functions  and  their  derivatives.  By  comparison  finding 
eigenvalues  by  the  first  approach  is  like  “shooting  in  the  dark". 

Another  advantage  of  the  second  approach  is  the  analytical  results  that 
can  be  obtained.  Consider  for  example  the  behavior  at  high 
frequencies.  The  eigenvalues  for  a  bounded  refractive  duct  were 
demonstrated  to  go  to  the  eigenvalues  for  an  unbounded  refractive  duct 
for  high  frequencies.  Similarly  the  eigenvalues  for  a  double  duct 
configuration  consisting  of  a  surface  duct  overlaying  a  refractive  duct 
were  demonstrated  to  go  to  the  composite  of  the  eigenvalues  for  a 
half-bounded  surface  duct  and  for  an  unbounded  refractive  duct. 

Another  example  of  an  analytical  result  is  the  development  of 
non-inteqral  mode  numbers.  With  the  use  of  non-integral  values  of  mode 
numbe’'  the  phase  integral  method  of  ray  theory  has  been  brought  into 


congruence  with  the  exact  solution  of  normal  mode  theory  for  the  case 
of  an  asymmetric  refractive  duct  without  boundaries.  The  ray  theory 
expressions  for  phase  and  group  velocity  are  identical  to  those  of  mode 
theory.  These  expressions  are  valid  for  all  frequencies.  They  are 
also  valid  for  sound  speed  profiles  in  which  attentuation  is  introduced 
by  means  of  complex  coefficients  in  the  profile  representation. 

A  third  example  of  an  analytical  result  is  the  presence  of  eigenvalues 

corresponding  to  pure  imaginary  frequencies.  Of  course  if  it  occurred 

to  one  to  do  so,  one  could  have  obtained  this  result  by  the  first 

2 

approach  by  letting  f  be  negative, 
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Table  1.  Comparison  of  the  exact  solution  of  Eq.  (27)  for 
moPe  1  with  three  algebraic  approximations. 


p 

Eq.  (27) 

Eq.  (35) 

Eq.  (37) 

Eq.  (42) 

1  .00 

1  .01879 

1.27599 

1 .01879 

1 .31784 

.95 

1  .07100 

1 .31096 

1  .07101 

1.36866 

.90 

1.12574 

1  .34645 

1.12576 

1  .42002 

.85 

1 .18295 

1  .38239 

1.18302 

1  .47179 

.80 

1  ,24254 

1  .41870 

1.24280 

1.52382 

.75 

1  .30436 

1  .45531 

1  .30510 

1.57596 

.70 

1  36823 

1  .49211 

1  .37001 

1  .62803 

.65 

1  .43393 

1.52902 

1  43774 

1  .67990 

.60 

1 .50120 

1.50592 

1.50870 

1  .73139 

.55 

1  .56975 

1  .60271 

1.58362 

1.78238 

.50 

1  .63928 

1.63928 

1  .66368 

1  .83273 

.45 

1.70949 

1  .67553 

1.75076 

1  .88237 

.40 

1.78009 

1 .71137 

1  .84778 

1  .93126 

.35 

1 .85083 

1  .74673 

1  .95931 

1  .97944 

.30 

1 .92148 

1  .78157 

2.09272 

2.02704 

.25 

1  ,99187 

1  .81591 

2.26050 

2.07434 

.20 

2.06189 

1.84984 

2.48560 

2.12181 

.15 

2.13148 

1  .88360 

2.81616 

2.17022 

.10 

2.20066 

1  .91763 

3.37727 

2.22079 

.05 

2.26948 

1  .95273 

4.66840 

2.27552 

.00 

2,33811 

1  .99037 

OD 

2.33811 

Table  2.  Counterpart  of  Table  1  for  mode  2. 


p 

Eq.  (27) 

Eq.  (35) 

Eq.  (37) 

Eq.  (42 

o 

o 

2.33811 

2.82734 

2.3381  1 

3.00431 

.95 

2.45472 

2.86793 

2.45471 

3.06841 

.90 

2.56867 

2.91018 

2.56935 

3.13205 

.85 

2.68156 

2.95431 

2.67978 

3.19507 

.80 

2.78946 

3.00058 

2.78346 

3.25731 

.75 

2.89295 

3.04932 

2.87767 

3.31861 

.70 

2.99201 

3.10094 

2.95963 

3.37880 

.65 

3.08688 

3.15599 

3.02668 

3.43771 

.60 

3  17792 

3.21517 

3.07646 

3.49517 

.55 

3.26554 

3.27944 

3.10712 

3.55102 

.50 

3.35012 

3.35012 

3.11748 

3.60511 

.45 

3.43199 

3.42908 

3.10712 

3.65732 

.40 

3,51143 

3.51911 

3.07646 

3.70759 

.35 

3.58870 

3.62448 

3.02668 

3.75592 

.30 

3.66402 

3.75214 

2.95963 

3.80239 

.25 

3.73761 

3.91428 

2.87767 

3.84729 

.20 

3.80968 

4.13424 

2.78346 

3.89111 

.15 

3.88048 

4.46324 

2.67978 

3.93473 

.10 

3.95025 

5.03981 

2.56935 

3.97967 

.05 

4.01930 

6.41385 

2.45471 

4.02870 

.00 

4.08795 

15.4774 

2.33811 

4.08795 

Table  3.  Roots  of  Eq.  (27)  and  relationship  to  the  non-integral 

mode  number  of  Eq.  (49). 


Mode  Modified 

Number  n 

1  .93643 

2  2.01734 

3  2.98458 

4  4.00790 


Upper  Duct 

X  Eq.  (49) 

1.509857  .976059 

3.189256  1.961474 

4.542342  2.984444 

5.713857  4.004702 


Lower  Duct 

X  Eq.  (49) 

1.234323  .943822 

2.775082  1.996165 

3.917366  3.009325 

4.879389  3.988301 
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Table  4.  Lower  order  roots  of  Eqs.  (65)  and  (68). 


Root 

Number 

Eq.  (65) 

Eq.  (68) 

1 

0 

0.4899060 

2 

1 .7647488 

1 .5621030 

3 

2.8082340 

2.3381074 

4 

3.6816163 

2.9624100 

5 

4.4606953 

3.5439131 
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Table  5.  Simplifications  arising  from  various 
conditions  between  the  parameters  of  Fig.  ll. 


Case  Condition  Variables 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 


General 


2  5 

3  5 

4  +  5 
2  +  3 
2+4 
3+4 


X.  y,  w.  P.  P^ 

X.  y.  w, 

X,  y,  w,  p 

X,  y,  w,  p.  p^ 

X.  w,  p,  p^ 

X.  w. 

X,  w.  p 
X.  w,  p 
X,  y,  w 
X.  y,  w,  p^ 

X.  y.  w,  p 

X,  W 
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FIG. 4.  Counterpart  of  the  curves  of  Fig.  3  for  non-integral  values  of 
mode  number 

FIG. 5.  Comparison  of  the  modified  phase  Integral  results  for  single 
ducts  with  the  normal  mode  result  for  a  double  duct. 

FIG. 6.  Schematic  of  the  bounded  refractive  duct. 

FIG. 7,  The  first  five  roots  of  Eq.  (55)  for  . 

FIG.0.  The  curve  segments  (of  Fig.  7),  which  lie  to  the  right  of  the 

s’ant  line,  represent  real  frequencies  Those  to  the  right  represent 

pure  Imaginary  frequencies. 

FIG. 9.  The  extension  of  the  first  root  of  Fig.  7  to  negative  values  of 

X  . 

FIG.  10.  The  first  four  roots  of  Eq  (55)  for  y=0. 
f^IG.  '1  Schematic  of  a  double  duct  profile 

►IS  '?  Schematics  of  bounded  ducts  lor  test'nq  the  ray  theor-e<,  o* 

5  and  6. 
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FIG.1.  Schematic  of  the  unbounded  refractive  duct. 
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FIG. 2.  The  eigenvalue,  x,  as  a  function  of  p  for  the  first  four 
modes . 
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FIG. 3.  Ttie  circles  reoresent  noneal  mode  peace  velocities  for  single 
onoounded  ducts.  The  curves  are  the  'ay-theory  phase-integral  results. 
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FIG. 4.  Counterpart  of  the  curves  of  Fig.  3  for  non-integral  values  of 
mode  number. 
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FIG. 5.  Comparison  of  the  modified  phase  integral  results  for  single 
ducts  with  the  normal  mode  result  for  a  double  duct. 
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G.6.  Schematic  of  the  bounded  refractive  duct 
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FIG. 7.  The  first  five  roots  of  Eq.  (55)  for  . 
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FIG. 8.  The  curve  segments  (of  F1g.  7),  which  lie  to  the  right  of  the 
slant  line,  represent  real  frequencies.  Those  to  the  right  represent 
pure  Imaginary  frequencies. 
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FIG. 9.  The  extension  of  the  first  root  of  F1g.  7  to  negative  values  of 

X  . 
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FIG.  10.  The  first  four  roots  of  Eq.  (55)  for  y=0. 
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